function C = slcov(X, w, vmean)
%SLCOV Compute the covariance matrix
%
% [ Syntax ]
%   - C = slcov(X) 
%   - C = slcov(X, w)
%   - C = slcov(X, w, vmean)
%   - C = slcov(X, w, 0)
%
% [ Arguments ]
%   - X:        the sample matrix
%   - w:        the weights of the samples
%   - vmean:    the pre-computed mean vector
%   - C:        the computed covariance matrix
%
% [ Description ]
%   - C = slcov(X) computes the covariance matrix for the samples in X.
%     the samples are stored as column vectors. 
%
%     If there are n samples in d-dimensional space, then X should be 
%     a d x n 2D matrix, where d is the vector dimension, and n is the 
%     number of samples. 
%
%   - C = slcov(X, w) computes the weighted covariance matrix for
%     the samples in X. 
%
%     As in slmean, w can be given in either of the following forms:
%       - a 1 x n row vector giving the sample weights. 
%       - a K x n matrix giving k groups of weights over all samples
%       - empty, simply meaning not weighting.
%
%   - C = slcov(X, w, vmean) computes the (weighted) covariance matrix
%     with the mean vector supplied. Thus in the function, vmean will be
%     used, instead of re-computing the mean vector. 
%
%   - C = slcov(X, w, 0) computes the (weighted) covariance matrix on 
%     the centered vectors. Since the vectors are treated as centered,
%     no mean vector would be computed, and X will not be shifted.
%
% [ Remarks ]
%   - M should be a 2D matrix (d x n), then w should be a 1 x n row vector,
%     vmean should be a d x 1 column vector or 0. Then C would be a 
%     d x d matrix.
%
%   - The function gives the maximum likelihood estimation of the
%     covariance matrix.
%
% [ History ]
%   - Created by Dahua Lin on Apr 22, 2006
%   - Modified by Dahua Lin on Sep 10, 2006
%       - replace sladd by sladdvec to increase efficiency
%       - replace slmul by slmulvec
%   - Modified by Dahua Lin on Jul 16, 2007
%       - base on new MATLAB function bsxfun
%   - Modified by Dahua Lin on Oct 18, 2007
%       - change the condition of NOT doing shifting to be
%         all(vmean == 0).
%   - Modified by Dahua Lin, on Dec 19, 2007
%       - add the functionality of computing covariance matrices for
%         multiple groups of weights.
%

%% parse and verify input arguments

error(nargchk(1, 3, nargin));

assert(isnumeric(X) && ndims(X) == 2, 'sltoolbox:slcov:invalidarg', ...
    'X should be a 2D numeric matrix.');
[d, n] = size(X);

if nargin < 2
    w =[];
    k = 1;
else
    if isempty(w)
        k = 1;
    else        
        assert(isnumeric(w) && ndims(w) == 2 && size(w,2) == n, ...
            'sltoolbox:slcov:invalidarg', ...
            'w should be a k x n numeric matrix if it is non-empty.');
        k = size(w, 1);
    end
end

if nargin < 3
    vmean = [];
else
    if ~isempty(vmean)        
        assert(isnumeric(vmean) && ...
            (isequal(vmean, 0) || isequal(size(vmean), [d, k])), ...
            'sltoolbox:slcov:invalidarg', ...
            'vmean should be either 0 or a d x k matrix if it is non-empty.');
    end
end

%% compute

% compute the mean vector

if isempty(vmean)
    vmean = slmean(X, w);
end

% compute covariances

if k == 1
    C = calccov(X, w, vmean);
else
    C = zeros(d, d, k);
    if isequal(vmean, 0)
        for i = 1 : k
            C(:,:,i) = calccov(X, w(i,:), 0);
        end
    else
        for i = 1 : k
            C(:,:,i) = calccov(X, w(i,:), vmean(:,i));
        end
    end
end

%% core computation

function C = calccov(X, w, vm)

% centralize the vectors
if ~slisallzeros(vm)
    Z = bsxfun(@minus, X, vm);
else
    Z = X;
end

% compute covariance
if isempty(w)
    C = ((1/size(X,2)) * Z) * Z';
else
    C = bsxfun(@times, Z, w / sum(w)) * Z';
end


